maxt <-read_csv("https://raw.githubusercontent.com/numbats/ida2024s2/master/data/vic_bushfire.csv") %>%pull(maxt)p1 <-ggplot() +geom_boxplot(aes(maxt)) +xlab(expression(degree*C)) +theme_light()p2 <-ggplot() +geom_density(aes(maxt)) +xlab(expression(degree*C)) +theme_light()patchwork::wrap_plots(p1, p2) + patchwork::plot_annotation(title ="Maximum daily temperature for bushfires in Victoria between 2015 and 2018")
Bivariate Relationships
When working with two numerical variables, a scatterplots is a natural choice to explore association between them.
It helps identify whether the variables are:
Positively associated: as one variable increases, the other also increases
Negatively associated: as one variable increases, the other decreases
No association: the data points are scattered randomly, showing no clear relationship
read_csv("https://raw.githubusercontent.com/numbats/ida2024s2/master/data/vic_bushfire.csv") %>%ggplot() +geom_point(aes(mint, maxt), alpha =0.6) +xlab(expression(Minimum~temperature~(degree*C))) +ylab(expression(Maximum~temperature~(degree*C))) +theme_light() +ggtitle("Daily maximum and minimum temperature for bushfires \nin Victoria between 2015 and 2018")
Correlation
Correlation is the linear association between two variables, ranging from \(-1\) to \(1\).
For two variables, \(X\) and \(Y\), the correlation coefficient \(r\) is calculated as follows:
where \(S_{xy}\) represents the biased sample covariance between \(X\) and \(Y\), and \(S_x\) and \(S_y\) denote the biased standard deviations of \(X\) and \(Y\), respectively.
Tip
In correlation, \(x\) and \(y\) are interchangeable, i.e. \(r_{xy} = r_{yx}\).
Correlation
Interpretation:
\(r = 1\): strong positive correlation
\(r = -1\): strong negative correlation
\(r = 0\): variables are not linearly associated, though this doesn’t imply they have no association!
Causation means that one variable directly influences another.
Here are some examples:
Smoking increases the risk of lung cancer.
Drinking alcohol impairs motor skills and cognitive function.
Remember: Correlation does not imply causation!
Both ice cream sales and drowning incidents rise in hot weather, but ice cream sales do not cause drownings. The real factor is the temperature.
As shoe size increases, so does a child’s reading ability. However, this relationship is actually due to age.
The number of pirates and global temperatures show a correlation, but this is a spurious relationship — two unrelated variables that happen to coincide!
What is Modeling?
Modeling is a method for understanding the relationship between different variables.
We typically represent the relationship of one response variable, \(Y\), in relation to other explanatory variables, \(X_1, X_2, \ldots, X_p\), using a mathematical function:
\[Y = f(X_1, X_2, ..., X_p).\]
Why we need modeling?
While we can observe the data, we may not have insights into the underlying processes, i.e. the joint distributions of the variables.
Therefore, we rely on data to make inferences about these processes and build models to explain them.
read_csv("https://raw.githubusercontent.com/numbats/ida2024s2/master/data/vic_bushfire.csv") %>%ggplot() +geom_point(aes(mint, maxt)) +geom_smooth(aes(mint, maxt), method ="lm", se =FALSE) +xlab(expression(Minimum~temperature~(degree*C))) +ylab(expression(Maximum~temperature~(degree*C))) +theme_light() +ggtitle("Daily maximum and minimum temperature for bushfires \nin Victoria between 2015 and 2018")
Simple Linear Regression
A simple linear regression (SLR) is a regression model that involves a single explanatory variable. The formula is expressed as follows:
\(x_i\) and \(y_i\): the \(i\)-th observation’s \(x\) and \(y\) value
\(\beta_0\) and \(\beta_1\): unknown constants that represent the intercept and slope, also known as coefficients or parameters.
\(\varepsilon\): the error term, which has specific properties and this is the stochastic component of the model, enabling us to make statistical inferences about the model and its parameters!
\[\hat{\beta}_0 = \bar{y} - \hat{\beta}_1\bar{x} \Rightarrow \text{intercept},\] where \(\bar{y} = \frac{1}{n}\sum_{i=1}^ny_i\) and \(\bar{x} = \frac{1}{n}\sum_{i=1}^nx_i\) are the sample mean.
For \(\text{mint} = 10\), \(\widehat{\text{maxt}} = 10 \times 1.04 + 12.02 = 22.42\). So, on average, we expect a day with a minimum temperature of 10 to have a maximum temperature of 22.42.
Warning: While we “expect” this to happen, there will be some variability, which is captured by the error term.
Predictions are generally valid for new values within the range of the original data used to fit the model. Extrapolation beyond this range is less reliable and carries greater risk.
Coefficient Interpretation
Note that, when \(\text{mint} = 0\), the predicted \(\widehat{\text{maxt}} = 0 \times 1.04 + 12.02 = 12.02\).
This value is equal to the intercept, indicating that the intercept represents the predicted average maximum temperature when the minimum temperature is zero.
Similarly, the slope \(\hat{\beta}_1 = 1.04\) can be interpreted as follows:
For each additional 1 degree increase in minimum temperature, the maximum temperature is expected to increase, on average, by 1.04 degrees.
Classical Normal Linear Regression Model
We have discussed the most basic regression model in statistics. Now, let’s consider its extension to include multiple explanatory variables:
where \(\varepsilon_i \sim N(0, \sigma^2)\) for \(i = 1,...n\).
Model assumptions:
Linearity: The relationship between explanatory variables and response variable is linear.
Independence: The errors are assumed to be independent of each other.
Homoscedasticity: The errors have constant variance across all levels of the explanatory variables.
Normality: The errors are assumed to be normally distributed.
Classical Normal Linear Regression Model
Coefficient Interpretation
In a CNLRM, the interpretation of the intercept term \(\hat{\beta}_0\) remains the same.
It represents the predicted average value of \(y\) when all \(x\) values are equal to zero.
For the other parameters \(\hat{\beta}_j\) where \(j = 1, \ldots, p\), the interpretation is as follows:
While holding all other variables constant, for each additional unit increase in \(x_j\), the variable \(y\) is expected to increase, on average, by \(\hat{\beta}_j\) units.
Let’s have a break!
Model Diagnostics
Model Diagnostics
All models are wrong; some are useful - George Box.
After fitting a linear model, it is crucial to assess whether the model assumptions are satisfied.
There is a real risk, that a model is imposing structure that is not really there.
It can result in inaccurate estimates of the coefficients and misleading interpretations of the relationship between \(x\) and \(y\).
Residual plots
The first thing to check is the residuals versus fitted values plot.
What constitutes a good residual plot?
Residuals should be randomly scattered around the horizontal line at zero.
There should be no discernible patterns in the residuals, such as curves or trends.
The spread of the residuals should remain consistent across all levels of the fitted values.
Residual versus Explanatory Variable
The second thing to check is the residuals versus each individual explanatory variable plot.
The principles for a good residual plot remain the same.
This type of plot offers clearer insights into the relationship between the residuals and each individual explanatory variable, which can be invaluable for refining and improving the model.
Q-Q plot
The third thing to check is the Q-Q plot.
To check a Q-Q plot, compare the distribution of your residuals to a theoretical normal distribution.
In a good Q-Q plot, the points should closely follow the 45-degree reference line.
Significant deviations from this line, such as curves or outliers, indicate that the residuals may not be normally distributed, which could suggest violations of model assumptions.
Lineup of Residual Plots
We can display the actual residual plot alongside a matrix of residual plots simulated under the assumption that the model is correct.
If the true residual plot is distinct and identifiable within the matrix, it provides evidence suggesting a violation of the model assumptions.
# remotes::install_github("autovi")mod <-lm(maxt ~ mint + se, data = bushfire)autovi::auto_vi(mod)$plot_lineup(alpha =0.1)
Model Selection
Beyond checking model assumptions, we often need to select between different models. For instance, choosing between maxt ~ rf + mint or maxt ~ mint.
There are several statistics to assess model fit:
\(R^2\)
adjusted \(R^2\)
AIC
BIC
Deviance
…
\(R^2\)
\[R^2 = 1-\frac{\sum_{i=1}^{n}(y_i-\hat{y})^2}{\sum_{i=1}^{n}(y_i-\bar{y})^2} = \frac{\text{Residual sum of square}}{\text{Total sum of square}} \in [0,1].\]
The coefficient of determination, \(R^2\), ranges from 0 to 1, with 1 indicating a perfect fit.
Adding more variables increases \(R^2\) but also adds model complexity.
Adjusted \(R^2\) penalizes for excessive variables to account for model complexity.
\(R^2\) is a common measure of linear model fit, indicating the percentage of variability in \(y\) explained by the model, with the remainder attributed to factors not included.
\(R^2\)
Which model is better?
AIC, BIC, and Deviance
AIC, BIC, and Deviance are goodness of fit measure to compare models.
AIC = Akaike Information Criterion (can be used to compare models. The smaller the value the better the model.)
Similarly BIC = Bayes Information Criterion indicates how well the model fits, best used to compare two models. Lower is better.
Deviance is the residual variation, how much variation in response that IS NOT explained by the model. The close to 0 the better, but it is not on a standard scale. In comparing two models if one has substantially lower deviance, then it is a better model.
Going beyond A Single Model
Gapminder
Hans Rosling was a Swedish doctor, academic and statistician, Professor of International Health at Karolinska Institute. Sadly he passed away in 2017.
He developed a keen interest in health and wealth across the globe, and the relationship with other factors like agriculture, education, energy.
